How to use skyscapeR
Fabio Silva
2018-12-13
skyscapeR is an open source R package for data reduction, visualization and analysis in skyscape archaeology, archaeoastronomy and cultural astronomy. It is intended to be a fully-fledged, transparent and peer-reviewed package offering a robust set of quantitative methods while retaining simplicity of use.
It includes functions to transform horizontal (Az/Alt) to equatorial (Dec/RA) coordinates, create or download horizon profiles, plot them and overlay them with visible paths of common celestial objects/events for prehistoric or historic periods, as well as functions to test statistical significance and estimate stellar visibility and seasonality. It also includes the ability to easily construct azimuth polar plots and declination curvigrams. Future versions will add more data reduction and likelihood-based model selection facilities, as well as an easy-to-use Graphical User Interface.
1. First Steps
Like all vignettes this is neither intended for those unexperienced in R nor a fully-fledged and detailed manual for skyscapeR. Neither will it be understandeable to those unfamiliar with the terminology and methodologies of cultural astronomy, archaeoastronomy or skyscape archaeology.
This document is an entry-point to the package’s main functions and sets out workflows of what you can do with it. More detail, as well as more functions, can be found in the package’s manual pages.
If you are new to R then I recommend the online free short course Code School’s TryR. For those who are already familiar with R, I apologise in advance for any patronising bits (such as the next two sections). With some basic R skills, guidance from this vignette and autonomous exploration of the package’s manual pages, anyone with prior training in skyscape archaeology can become a master skyscapeR.
Installation
The package requires that the latest version of R is installed first. See the R Project website for details. Also suggested is the installation of RStudio. With R installed, you can install the latest skyscapeR release directly from CRAN by doing:
install.packages('skyscapeR')
This will install the package itself along with any dependencies. Upon successful completion you should see a line saying * DONE (skyscapeR). If this doesn’t happen then there should be an error message explaining what went wrong.
If you rather install the latest development version (often unstable or untested) you can do so directly from the GitHub repository. Firstly, ensure you have package devtools installed, by running:
install.packages('devtools')
Then type the following to download and installthe development version of skyscapeR:
devtools::install_github('f-silva-archaeo/skyscapeR')
At the end, you should see the same successful completion line.
Initialization
Every time you want to use the package it needs to be loaded. For this you need to type:
## Loading required package: swephR
## Warning: replacing previous import 'httr::config' by 'plotly::config' when
## loading 'skyscapeR'
## Warning: replacing previous import 'graphics::layout' by 'plotly::layout'
## when loading 'skyscapeR'
## Warning: replacing previous import 'plotly::filter' by 'stats::filter' when
## loading 'skyscapeR'
Datasets and help
The current version of skyscapeR comes with a couple of datasets of measurements that can be used for learning and testing the package. These can be loaded by typing:
data(RugglesRSC)
or
data(RugglesCKR)
If you want to know more about these datasets then use the helpful ? command which opens up the help page related to whatever you are asking help for. In this case, if you want to know more about the RugglesRSC dataset you can type in the console ?RugglesRSC. This opens the manual page for the dataset, or function. In RStudio this opens on the bottom-right pane by default.
2. Hello Curvigram!
Let’s start by creating a histogram with skyscapeR. The function for this is called histogram(). Go ahead and do ?curvigram right now to learn more about it, including an example at the bottom of the manual page, which we will now do:
data(RugglesRSC)
hist <- histogram(RugglesRSC$Dec, 2)
This creates a histogram based on the declination data in the RugglesRSC dataset, and using an uncertainty of 2º for all measurements. You can visualize it by typing:
Be sure to check ?histogram to see what other options are available.
A curved histogram, also known as curvigram, can be created from the same data by using curvigram() in a similar manner to histogram(). This is a more robust form of data visualization that can use different measurement uncertainties for the different measurements. In the next example a set of random uncertainties (variable unc) are used to construct a curvigram:
data(RugglesRSC)
unc <- runif(length(RugglesRSC$Dec),1,10) # random values between 1 and 10
curv <- curvigram(RugglesRSC$Dec, unc)
plot(curv)
Check ?curvigram to see what other options are available.
Adding celestial objects to a curvigram or histogram
Let’s add some celestial targets to the curvigram on order to compare them with those frequency peaks. We need to create a skyscapeR.object first. This is done with sky.objects():
lunar <- sky.objects('moon', epoch=-2000, col='red', lty=2)
This creates an object that includes all standard lunar targets (currently only the lunar extremes), set for the year 1999 BCE (since there is no year zero). These will be drawn in red colour with line type (lty) two, which is to say as dashed lines (check out other types and options in ?par). Then redo the curvigram with this object:
You can now see the southern major lunar extreme (sMjLX) and the southern minor lunar extreme (smnLX) declinations for the year 1999 BCE on the same plot.
3. Handling and Visualizing Fieldwork Data
skyscapeR has some built-in function to automatically process fieldwork data. Although the data reduction process can be done manually using other functions (such as az2dec(), mag.dec() or sunAz()) the entire process has been automated by using the reduct. functions.
Data reduction of compass measurements
Imagine that you’ve just returned from a day of fieldwork surveying a set of archaeological structures for their orientation with your compass and clinometer. You can input your measurements into R, along with the georeferences and measurement data, and use reduct.compass() to process it:
georef <- rbind( c(35.1, -7.1), # GPS data
c(35.1, -7),
c(35.2, -7.1),
c(35.1, -7.3) )
azimuths <- c(93, 108, 105, 98) # Compass measurements
altitudes <- c(2, 1.5, 0.5, 1) # Clinometer measurements
data <- reduct.compass(loc=georef, mag.az=azimuths, date="2017/06/13", alt=altitudes)
## Altitude values found. Calculating declination...
## Latitude Longitude Magnetic.Azimuth Date Mag.Dec True.Azimuth
## 1 35.1 -7.1 93 2017/06/13 -1.704983 91.295
## 2 35.1 -7.0 108 2017/06/13 -1.674106 106.326
## 3 35.2 -7.1 105 2017/06/13 -1.705037 103.295
## 4 35.1 -7.3 98 2017/06/13 -1.766996 96.233
## Altitude Declination
## 1 2.0 0.091
## 2 1.5 -12.407
## 3 0.5 -10.552
## 4 1.0 -4.518
This has: (1) automatically retrieved the value of magnetic declination for your locations and date; (2) corrected azimuths to true; and (3) calculated the corresponding declination. The resulting data.frame object can be easily exported into a file for safe keeping using the standard R utils (such as write.csv(), for example).
If the altitude component is not included then the function will only calculate the true azimuth:
data <- reduct.compass(loc=georef, mag.az=azimuths, date="2017/06/13")
## No altitude values or horizon profile found. Declination values were not calculated.
## Latitude Longitude Magnetic.Azimuth Date Mag.Dec True.Azimuth
## 1 35.1 -7.1 93 2017/06/13 -1.704983 91.295
## 2 35.1 -7.0 108 2017/06/13 -1.674106 106.326
## 3 35.2 -7.1 105 2017/06/13 -1.705037 103.295
## 4 35.1 -7.3 98 2017/06/13 -1.766996 96.233
In such cases, if a horizon profile is given, the altitude can be automatically retrieved from it. See section on Dealing with Horizons below for examples.
Data reduction of theodolite measurements
On the other hand, if fieldwork consisted on theodolite/total station measurements using the sun-sight technique, the azimuths still need to be processed. For such precise measurements we recommend using function ten() to convert from deg-arcmin-arcsec to decimal point degree as follows:
georef <- c( ten(35,50,37.8), # GPS data
ten(14,34,6.4) )
az <- c( ten(298,24,10), # Theodolite H measurements
ten(302,20,40) )
alt <- c( ten(1,32,10), # Theodolite V measurements
ten(0,2,27) )
az.sun <- ten(327,29,50) # The azimuth of the sun as measured at time
limb <- "right" # Which limb of the sun was targeted
date <- "2016/02/20"
time <- "11:07:17" # Time the sun measurement was taken
timezone <- "Europe/Malta" # Timezone corresponding to time above
data <- reduct.theodolite(loc=georef, az, date, time, timezone, az.sun, limb, alt)
## Altitude values found. Calculating declination...
## Latitude Longitude Uncorrected.Azimuth Date.Time Sun.Az
## 1 35.84383 14.56844 298.4028 2016-02-20 11:07:17 185.9292
## 2 35.84383 14.56844 302.3444 2016-02-20 11:07:17 185.9292
## True.Azimuth Altitude Declination
## 1 156.8347 1.53611111 -46.82722
## 2 160.7764 0.04083333 -49.90692
Similarly to reduct.compass() if a horizon profile is given, the altitude can be automatically retrieved from it.
Plotting azimuths
In general, the use of azimuths for analysis and visualization in skyscapeR is deprecated since azimuths are location-specific. It is preferable to convert all measurements to declinations and work with equatorial coordinates (see section below). However skyscapeR does include a much requested function to create a polar plot of azimuth values. The function is plotAz():
az <- rnorm(30, 85, 20) # This creates 30 random azimuths
plotAz(az)
You can use the same skyscapeR.object for this plot, but then you need to specify a single location, since azimuths are location-specific. At the moment the horizon altitude is assumed to be 0º and flat as well.
sunandmoon <- sky.objects(c('sun','moon'), epoch=-4000, col=c('black','red'), lty=c(2,3), lwd=2)
plotAz(az, obj=sunandmoon, loc=c(52,0))
Converting azimuth to declination
If you are looking for a mere horizontal to equatorial coordinate conversion, without all the extra automation that the reduct. functions provide you can use the az2dec() function:
dec <- az2dec(az, loc= c(35,-7), alt=0)
dec
## [1] 3.7546636 -1.8415784 31.0691609 -12.0497851 -1.1302515
## [6] -10.8423548 7.3101096 -7.2664169 0.1768612 19.1757410
## [11] 3.3432657 -27.8503957 12.5812751 -21.2983433 -1.7473546
## [16] -11.4197380 -19.3548093 19.0381906 10.9146679 5.3136124
## [21] 23.1304054 4.9296942 -3.0214476 2.7376682 16.2262860
## [26] 4.1979081 10.9197771 22.7497600 5.7338989 12.0847213
4. Handling Horizon Profiles
One of skyscapeR’s main missions is to handle horizon profile data. This is done through the creation of a skyscapeR.hor object which can then be plugged in to several functions of this package not only for automation, but also for the purposes of visualizing the visible paths of celestial objects.
Creating or downloading profiles
If you have azimuth and altitude data a profile can be constructed as in the following example:
az <- c(0,90,180,270,360)
alt <- c(0,5,5,0,0)
georef = c(40.1, -8)
hor <- createHor(az, alt, loc=georef, name= 'Horizon Profile 1')
This can be visualized by simply typing:
Not the prettiest horizon you’ve ever seen, but that’s what can be interpolated from the five datapoints given…
A prettier horizon profile can be downlaoded from the HeyWhatsThat (HWT) website. These profiles are based on the SRTM digital elevation models and, therefore are not always trustworthy (especially so if the horizon is near). Nevertheless they can be a great help for situations when a horizon altitude is impossible to measure on site. To do this, first create a horizon profile at the HWT website, then save the 8-digit ID code that is as the end of the permanent link given by HWT. For example, if the link is https://www.heywhatsthat.com/?view=NML6GMSX, save the bit after view= and use function downloadHWT():
hor <- downloadHWT('NML6GMSX')
plot(hor)
Exporting profiles to Stellarium
One can then export these profiles into a format that Stellarium recognises:
exportHor(hor, name='Test', description='Test horizon export to Stellarium')
This creates a zip-file ready to be imported into Stellarium’s landscape feature.
Automating horizon altitude retrieval
With a horizon profile set up (or even many), you can retrieve the horizon altitude for any azimuth value by doing:
## alt alt.unc
## 1 1.79 0.02
## alt alt.unc
## 1 2.73 0.02
This automation can also be used in the data reduction functions by using a skyscapeR.horizon object rather than simple georeferences:
data <- reduct.compass(loc=hor, mag.az=az, date="2017/06/13")
## Horizon profile found. Obtaining altitude values and calculating declination...
## Latitude Longitude Magnetic.Azimuth Date Mag.Dec True.Azimuth
## 1 40.44385 -7.938178 0 2017/06/13 -2.118487 -2.118
## 2 40.44385 -7.938178 90 2017/06/13 -2.118487 87.882
## 3 40.44385 -7.938178 180 2017/06/13 -2.118487 177.882
## 4 40.44385 -7.938178 270 2017/06/13 -2.118487 267.882
## 5 40.44385 -7.938178 360 2017/06/13 -2.118487 357.882
## Altitude.alt Altitude.alt.unc Declination
## 1 1.29 1.37 50.799
## 2 1.99 0.02 2.902
## 3 1.20 0.02 -48.311
## 4 0.81 1.50 -1.086
## 5 1.29 1.37 50.799
Visualizing celestial object paths
plot() can be used to display the visible paths of celestial objects chosen with sky.objects(). To use the one we created before for the curvigram just type:
Or create a new one, in this case using an epoch range for cases where we might have uncertainty in the age of the site, and including both solar and a stellar target:
aux <- sky.objects(names=c('sun', 'Aldebaran'), epoch=c(-4300,-3700), col=c('blue', 'red'))
plot(hor, obj=aux)
5. Stars
skyscapeR includes data on the brightest 100 stars in the night sky. The data table is accessible via data(stars). To pick a particular star, maybe to see what declination it has, or had, you can do as follows:
## $name
## [1] "Sirius"
##
## $constellation
## [1] "alCMa"
##
## $app.mag
## [1] -1.46
##
## $ra
## [1] 101.2842
##
## $dec
## [1] -16.71739
##
## $epoch
## [1] 2000
##
## attr(,"class")
## [1] "skyscapeR.star"
This shows the information for Sirius which is now loaded into object ss. Note the epoch the data is output, J2000, by default. To get coordinates for other epochs just do:
ss <- star('Sirius', year=-4000)
ss$dec
## [1] -26.40989
Star phases and seasonality
To estimate the phases, events and/or seasonality of stars, the function star.phases() can be used. It works as follows, for the location of Cairo, an epoch of -3000, and a horizon altitude of 2º:
# ncores forced to 2 for production of this document
sp <- star.phases('Sirius', -3000, loc=c(30.0,31.2), alt.hor=2, ncores=2)
## Running calculations on 2 processing cores. This may take a while...
One can then check the star’s phase, event date-range and seasons by typing:
## [1] "Arising and Lying Hidden"
## [,1]
## Achronycal Setting "March 8 - March 30"
## Heliacal Rising "June 5 - June 29"
## [,1]
## Rising "June 5 - September 23"
## Rising and Setting "September 24 - December 4"
## Setting "December 5 - March 30"
## Arising and Lying Hidden "March 31 - June 4"
The time period for events is a range of dates where, given the parameters, it is possible to observe the star and recognise the observation as the said event (note in particular the meaning of parameter alt.rs in the manual page for star.phases()). These can be visualized using:
Significance Testing for Curvigrams
You can test the statistical significance of an empirical curvigram by comparing it with the expectattion of a given null hypothesis. This technique can be used to output a p-value which is an established measure of significance.
To demonstrate significance testing we will again use the Recumbent Stone Circle data. To do this one uses sigTest() to run the significance testing routine. This can take a while depending on your machine’s resources. If it takes too long, try lowering the nsims parameter (though this brings a cost of resolution, see the manual page for this function).
data <- data.frame(True.Azimuth = RugglesRSC$CL_Rec_C)
data$Azimuth.Uncertainty <- rep(5, length(data$True.Azimuth))
data$Latitude <- RugglesRSC$Latitude
data$Altitude <- RugglesRSC$CL_Mean_Alt
# ncores forced to 2 and nsims to 100 for production of this document
sg <- sigTest(data, ncores=2, nsims=100)
## Creating Empirical Curvigram...Done.
## Loading required package: foreach
## Running 100 simulations on 2 processing cores. This may take a while...Done.
## Performing a 2-tailed test at the 95% significance level.
One can then plot the curvigram again, but now with the results of the significance testing displayed. This adds the expectation around the null hypothesis (the grey shaded area) as well as the estimated global p-value.
One can highlight the regions of significance, i.e. those that are outside the region of significance by doing:
If one is not interested in ploting the results of the significance testing, but simply getting the values then you can retrieve them from the output of sigTest():
##
## *** Results of Significance Test ***
##
## 2-tailed test at 95% confidence, based on 100 simulations.
## global p-value: 0.0099 (**)
## local p-values:
## + dec range [-32.97, -27.51] :: p-value: < 0.01 (**)
## + dec range [-27.47, -26.74] :: p-value: < 0.01 (**)